function [spikes Gampa]=izhikevich_fs(stim_conductance_exc, stim_conductance_inh)

pyramid=true;
if pyramid
    C=104 % *pF
    gL=4.3 %*nS
    taum=C/gL
    EL=-65 %*mV # Same as changing I
    VT=-52 %*mV
    DeltaT=.8 %*mV
    Vcut=VT+5*DeltaT
    b=65 %*nA
    Vr=-53 %*mV
    tauw=88
    a=-.8
else
    % FS cell 
    C=59 % *pF
    gL=2.9 %*nS
    taum=C/gL
    EL=-62 %*mV # Same as changing I
    VT=-42 %*mV
    %VT=-32 %*mV
    DeltaT=3 %*mV
    Vcut=VT+5*DeltaT
    b=61 %*nA
    Vr=-54 %*mV
    tauw=16
    a=1.8

end

T=1000;

dt_ms=1;


spikes = 0;
n=round(T/dt_ms);

v=zeros(1,n);
u=zeros(1,n);
v(1)=EL;
u(1)=0;


stp_tau = exp(-dt_ms/150);
stp = 0*stp_tau.^(0:n-1)+1;

si= stim_conductance_inh .*[zeros(1,5) ones(1,n-5)].*stp;
se= stim_conductance_exc .*[zeros(1,5) ones(1,n-5)].*stp;

% model conductances.
tau_ampa  = exp(-dt_ms/5);
tau_nmda  = exp(-dt_ms/150);
tau_gabaa = exp(-dt_ms/6);
tau_gabab = exp(-dt_ms/150);

E_AMPA	= 0.0;
E_NMDA	= 0.0;
E_GABAa	= -70;
E_GABAb = -90;
Ggabaa = 0;Ggabab = 0; Gampa = 0;Gnmda = 0;

Ggabaa=zeros(1,n);
Ggabab=zeros(1,n);
Gampa=zeros(1,n);
Gnmda=zeros(1,n);
I=zeros(1,n);

gain_ampa=1;
gain_nmda=.1;
gain_gabaa=1;
gain_gabab=.1;
% end model conductances.


for i = 2:n
    
    % conductances from synapses.
    Ggabaa(i) = Ggabaa(i-1)*tau_gabaa + (1-tau_gabaa)*gain_gabaa*si(i);
    Ggabab(i) = Ggabab(i-1)*tau_gabab + (1-tau_gabab)*gain_gabab*si(i);
    
    Gampa(i) = Gampa(i-1)*tau_ampa + (1-tau_ampa)*gain_ampa*se(i);
    Gnmda(i) = Gnmda(i-1)*tau_nmda + (1-tau_nmda)*gain_nmda*se(i);

    xx = (-80-v(i))/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    
     
    %g= ( Gampa(i+1)       +NMDAgate.*Gnmda(i+1)        + Ggabaa(i+1)         + Ggabab(i+1));
	%E= ( Gampa(i+1)*E_AMPA+NMDAgate.*Gnmda(i+1)*E_NMDA + Ggabaa(i+1)*E_GABAa + Ggabab(i+1)*E_GABAb);
    %I(i) = v(i).*g-E;
    I(i) = (v(i-1)-E_AMPA)*Gampa(i) + (v(i-1)-E_NMDA)*NMDAgate*Gnmda(i)+ ...
        (v(i-1)-E_GABAa)*Ggabaa(i)+(v(i-1)-E_GABAb)*Ggabab(i);
    
    
    dvm=(gL*(EL-v(i-1))+gL*DeltaT*exp((v(i-1)-VT)/DeltaT)-I(i)-u(i-1))/C; % volt
    v(i) = v(i-1) + dt_ms*dvm;
    dw=(a*(v(i-1)-EL)-u(i-1))/tauw; % amp
    u(i)=u(i-1) + dt_ms*dw;
    
    if v(i) < -90
        v(i)=-90;
    end
    
    if( v(i) > Vcut )
        v(i-1)=0;
        v(i) = Vr;
        u(i) = u(i)+b;
        spikes = spikes + 1;
    end

   
end
figure(3)
subplot(3,1,1)
plot(dt_ms*(1:n),[v' u']);
subplot(3,1,2)
plot(dt_ms*(1:n),[Gampa' Gnmda' Ggabaa' Ggabab']);
title('conductances');
subplot(3,1,3);
plot(dt_ms*(1:n),I);
title('Isyn (pA)');
drawnow
